Experimental Quantum Computing to Solve Systems of Linear Equations 
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Solving linear systems of equations is ubiquitous in all areas of science and engineering. With 
rapidly growing data sets, such a task can be intractable for classical computers, as the best known 
classical algorithms require a time proportional to the number of variables A''. A recently proposed 
quantum algorithm shows that quantum computers could solve linear systems in a time scale of 
order log (A''), giving an exponential speedup over classical computers. Here we realize the simplest 
instance of this algorithm, solving 2x2 linear equations for various input vectors on a quantum 
computer. We use four quantum bits and four controlled logic gates to implement every subroutine 
required, demonstrating the working principle of this algorithm. 



ng, China 



Oh. 



> 
O 

O 

cn 



X 



The problem of solving a system of linear equations 
plays a central role in diverse fields such as signal pro- 
cessing, economics, computer science, and physics. Such 
systems often involve tera or even petabytes of data, and 
thus the number of variables N, is exceedingly large. 
However, the best known algorithms for solving a system 
of N linear equations on classical computers requires a 
time complexity on the order of iV, posing a formidable 
challenge. 

Harnessing the superposition principle of quantum me- 
chanics, quantum computers [1. 2] promise to provide 
exponential speedup over their classical counterparts for 
certain tasks. Notable examples include quantum sim- 
ulation [3, 4] and Shor's quantum factoring algorithm 
[5], which have driven the field of quantum information 
over the past two decades as well as generating signifi- 
cant interest in quantum technologies that have enabled 
experimental demonstrations of the quantum algorithms 
in different physical systems [6-11]. 

Recently, Harrow, Hassidim and Lloyd [12] proposed 
another powerful application of quantum computing for 
the very practical problem of solving systems of linear 
equations. They showed that a quantum computer can 
solve a system of linear equations exponentially faster 
than a classical computer in situations that we are only 
interested in expectation values of an operator associated 
with the solution rather than the full solution. A quan- 
tum algorithm has been designed such that the value of 
this property may be estimated to any fixed desired ac- 
curacy within O (log (TV)) time, making it one of the most 
promising applications of quantum computers. 

In this article, we report an experimental demonstra- 
tion of the simplest meaningful instance of this algorithm, 
that is, solving 2x2 linear equations for various input 
vectors. The quantum circuit is optimized and compiled 



into a linear optical network with four photonic quantum 
bits (qubits) and four controlled logic gates, which is used 
to coherently implement every subroutine for this algo- 
rithm. For various input vectors, the quantum computer 
gives solutions for the linear equations with reasonably 
high precision, ranging from fidelities of 0.825 to 0.993. 



THEORY 

The problem of solving linear equations can be summa- 
rized as follows: We aim to solve Ax = b for x, when given 
a. N X N Hermitian matrix A and a vector b. To adapt 
this problem to quantum processing, x and b are scaled 
to unit length (i.e., ||x|| = ||6|| = 1). Thus, a vector b 
can be represented by a quantum state \b) = bi \i) on 
0(log(A^)) qubits where |i) denotes computational basis. 
The desired solution x can then be encoded within the 
quantum state as 



\x) = cA-^ \b) 



\b) 



(1) 



The quantum algorithm devised in ref. [12] was designed 
to synthesize \x) (see Fig. lA). The quantum algorithm 
involves three subsystems: a single ancilla qubit initial- 
ized in |0), a register of n qubits of working memory 
initialized in |0)®" and an input state initialized in \b). 
The input state |5) can be expanded in the basis of \uj) 

1^) = f^j Wj)^ where \uj) is eigenstate of A, and 

= {uj\b). Execution of the algorithm can be decom- 
posed into three subroutines: (1) phase estimation, (2) 
controlled rotation and (3) inverse phase estimation. 

Step (1) is used to determine the eigenvalues of A, 
which we denote by Xj. Phase estimation is essentially 
a controlled unitary with a change of basis that maps 
the eigenvalues onto the working memory [13, 14]. The 
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FIG. 1: Quantum circuits for solving systems of linear equations. A. Outline of the original quantum algorithm proposed in [12]. 
The light gray blocks represent three basic subroutines of the algorithm. U — Yl'kZi |fe)(fc| (g)e''*'°'°''"^, where T — 2^ with t being 
the number of registers. H is a Hadamard gate. FT and FT"^ are the Fourier transformation and the inverse Fourier transforma- 
tion [1], respectively. The controlled rotation R evolves the system into the state /3i|wj)|Aj)-y((l — C^/Aj)|0) + (C/Aj)|l)). 
The inverse phase estimation subroutine restores the register to |0)®". Finally, the vector x, the solution of the system of linear 
equations, can be obtained by conditioning on the measurement outcome of |1) in ancilla qubit. B. The optimized circuit with 
four qubits and four entangling gates (see main text for details). The two swap gates are canceled out. 



phase estimation protocol is applied to the input, using 
the working memory as control, to give 



N 



(2) 



Here |Aj) represents the binary representation of Aj, 
stored to a precision of n bits. 

In step (2), one needs to extract the eigenvalues of 
i.e. A""'^ from |Ai). This is realized through an additional 
ancillary qubit initialized in the state |0). Application of 
an appropriate controlled rotation i?(A~^) on this qubit 
(see Fig. lA) transforms the system to 




i-fTlo> + f|i»- 

\j Aj 



(3) 



The final step involves applying the gate sequence of step 
(1) in reverse. This disentangles the register, which is 
reset to |0)®". Therefore we end up with 




i-fTlo> + ^|i». 



A, 



(4) 



Measurement of the ancillary qubit with an outcome of 
|1) will result in an output state C{/3j / Xj)\uj) which 
is proportional to our expected result state \x). 



EXPERIMENT 

Here we demonstrate a proof-of-principle experiment 
of this algorithm: solving systems of 2 x 2 linear equa- 
tions. We choose the matrix A to be 



A 



1.5 0.5 
0.5 1.5 



(5) 



The eigenvalues of this matrix are 1 and 2, and can there- 
fore be encoded with two qubits in the register. We 
choose the following values for input vector |6) 



\bi) 



1 



,\b2) 



1 

71 



1 

-1 



(6) 



This allows us to optimize the circuit requiring four 
qubits and four entangling gates as shown in Fig. IB. 
The phase estimation subroutine of the circuit can be 
compiled into two controUed-NOT (CNOT) gates, a swap 
gate, and three single qubit rotation gates. Following the 
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FIG. 2: Experimental setup. There are four key modules in the optical setup. (1) Qubit initialization: Ultraviolet laser pulses 
with a central wavelength of 394 nm, pulse duration of 120 fs and a repetition rate of 76 MHz pass through two /3-barium borate 
(BBO) crystals to produce two photon pairs. The four single photons are spatially separated by PBSi and PBS2 and initialized 
using HWPs to be in the state \H)i with i denoting their spatial modes. Photon 1 is used as the input vector qubit and photon 
4 is used as the ancilla. Photons 2 and 3 are used as the register qubits Rl and R2, respectively. (2). Phase estimation: The 
input qubit \b) is mixed with the two register qubits on PBS3 and PBS4 to simulate the CNOT gates in Fig. IB. (3) R(X~^) 
rotation: An entanglement-based implementation of the two controlled unitary gates (see main text for details). (4) Inverse 
phase estimation: This is realized semiclassically by using measurement and classical feed forward. To achieve good spatial and 
temporal overlap, all photons are spectrally filtered (Afwhw = 3.2 nm) and detected by fiber-coupled single-photon detectors 
(Dl, • ■ • , D4). The coincidence events are registered by a programmable multichannel coincidence unit. 



circuit design of Ref. [15], the R(A~"'^) rotation subrou- 
tine is implemented in two steps: finding the reciprocal 
|l/Aj) from eigenvalue \Xj) stored in registers, which in 
our case can be realized by a swap gate, and controlled 
unitary gates H{6), where 

- 1, sin(20) -cos{2e) )■ 

Finally, the subroutine of the inverse phase estimation is 
realized using a semiclassical version that employs single- 
qubit rotations conditioned on measurement outcomes 
[16]. 

To implement the quantum circuit shown in Fig. IB, 
we prepare four single photons from spontaneous para- 
metric down-conversion [17] as the input qubits (Fig. 2). 
The horizontal (H) and vertical (V) polarizations of the 
single photons are used to encode the logic qubits |0) and 
|1), respectively. The experimental challenge of imple- 
menting the circuit in Fig. IB lies in the four entangling 
gates between the single photonic qubits. 

In the phase estimation subroutine, noting that the 
target qubits of the CNOT gates are fixed, their im- 
plementations can be simplified using combinations of 
a polarization beam splitter (PBS) and a half-wave plate 
(HWP), through which an arbitrary control qubit a\H) + 
I3\V) and the target qubit \H) evolve into a\H)\H) + 
P\V)\V) which is the desired output of CNOT operations 



[7, 18]. The R(A^^) rotation subroutine involves two con- 
secutive controlled unitary gates, H(7r/8) and H(7r/16). 
Instead of decomposing it into multiple CNOT gates 
[1, 18], we adopt a more efficient, entanglement-based 
construction method [19]. The ancilUa qubit is first en- 
tangled with the register qubits by mixing on PBS5, and 
then passed through a polarization-dependent Sagnac- 
like interferometer where the desired controlled unitary 
operations are applied (see SOM text for more details). 
Finally, the ancillary qubit is measured, and when an 
outcome state |1) is obtained, the algorithm is announced 
successful. 



RESULTS 

Before running the algorithm, we first characterized 
the performance of the optical quantum circuit. The two 
registers, ancilla and input qubits (I63)) are initialized in 
the \H)a \H)ri \H)r2 «) \H)b state. Theoretically 
these four qubits will evolve into a maximally-entangled 
four-qubit Greenberger-Horne-Zeilinger (GHZ) state dur- 
ing the subroutine of phase estimation and R{X~^) rota- 
tion. After the four photons pass through PBS3, PBS4 
and PBS5, we observed the Hong-Ou-Mandel type in- 
terference among the four photons. We measured the 
fidelity - defined as the overlap of the experimentally 
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FIG. 3: Experimental results. Three different input vectors are chosen: A |6i) = {\H) + \V))/V2, B I&2) = {\H) - \V))/V2, 
and C I63) = \H). The quantum algorithm is run to determine the the expectation value (a;|M|a;), where M is some operator. 
For each input state \b), the theoretically predicted (gray bar) and experimentally measured (red bar) expectation values of 
the observables of the Pauli matrices Z, X, and Y are presented. The output states are measured with a fidelity of 0.993(3), 
0.825(13), and 0.836(16) for \b2), and Ifea), respectively. The error bars denote one standard deviation, deduced from 

propagated Poissonian counting statistics of the raw detection events. 



produced state with the ideal one [20] - of the generated 
four-photon GHZ state [21] (see SOM Fig. SI). The mea- 
surements yield a state fidelity of 0.65(1), which exceeds 
the threshold of 50% [22] by 15 standard deviations. This 
confirms the presence of genuine entanglement [23] cre- 
ated during the quantum computation. 

We have implemented the algorithm for various input 
vectors |6) which are varied by tuning the HWP in front 
of PBS3. In accordance with Fig. IB, the two registers 
should be projected to the state |0), the ancilla qubit to 
state |1). The output \x) is measured in some desired 
observable. In the experiment, each run of the algorithm 
is finished by a fourfold coincidence measurement where 
all four detectors fires simultaneously. 

We characterize the output by measuring the expec- 
tation values of the Pauli observables Z, X, and Y for 
each input state \b). Figure 3 shows both the ideal (gray 
bar) and experimentally obtained (red bar) expectation 
values for each observable. To quantify the algorith- 
mic performance, we compute the output state fidelity 
F = (x\px\x), where \x) is the ideal state and Px is the ex- 
perimentally reconstructed density matrix of the output 
state from the expectation values of the Pauli matrices 
(see SOM Fig. S2). Compared with ideal outcomes, the 
output states have fidelities of 0.993(3) for |&i), 0.825(13) 
for [62), and 0.836(16) for [63), respectively. 

The difference in the performance for the three inputs 
is linked to the specific optical setup used in the exper- 
iment. The fidelity imperfections for [62) and I&3) are 
caused by high-order photon emission events. However, 
in the case for high-order photon emissions do not 
give a negative contribution, giving rise to a near-ideal 



algorithm performance. 



DISCUSSION 

In summary, we have presented a proof-of-principle 
demonstration of the quantum algorithm for solving sys- 
tems of linear equations in a small-scale quantum com- 
puter involving four qubits and four entangling gates. We 
have implemented every subroutine at the heart of the 
algorithm and characterized the circuit and algorithmic 
performances by the quantum state fidelities. The tech- 
nique of coherently controlling multiple qubits and exe- 
cuting complex, multiple-gate quantum circuits presents 
an advance on linear optics quantum computation [21, 25] 
and allows to test other similar quantum algorithms such 
as solving nonlinear differential equations [26, 27] and 
data fitting [28]. 

In principle, efficient quantum computation can be 
achieved using single-photon sources, linear optics, and 
single-photon detectors [24, 25, 29]. The current ex- 
periment, however, is still limited by a probabilistic 
single-photon source and inefficient detectors. It can 
be expected that with ongoing progress on determin- 
istic single-photon sources [30], high-efficiency (>93%) 
single-photon detectors [31], and on-chip integration [9], 
a larger-scale quantum circuit for solving more complex 
linear equations can be implemented in the future. 
note: During the stage of manuscript preparation, we 
became aware of related work from Walther's and Du's 
group [32]. 
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Experimental details 

In the controlled unitary gate, the control qubit turns 
on or off the unitary operation acting on target qubit. If 
the controlled photon is in \H), nothing happens to the 
target photons, whereas if the control photons is \ V), the 
target photon undergoes the H{9) (see eqn. 7) rotation. 
Although it is possible to decompose the controlled uni- 
tary gates into a series of CNOT gates and single-qubit 
unitary gates [1], it can be experimentally inefficient and 
expensive. Here, we take advantage of the recently pro- 
posed schemes for simplifying the construction of con- 
trolled unitary gates [19]. 

As shown in Fig. 2, first, the ancillla qubit (photon 
4) is entangled with the register qubit R2 (photon 3) by 
mixing them on PBS5. The register qubit R2 is after- 
wards entangled with the input state \h) = fi\H) + 
and the other register qubit Ri. Thus, after passing 
through PBS3, PBS4, PBS5, and single qubit gate X and 
H, the quantum state of the four qubits can be written as: 

ili+iy)\+)t\H)m\V)R2\H)A+{li-iy)\-)b\V)m\H)R2\V)A 
(normalization omitted). 

The ancillary qubit then passes through PBSg, where 
two spatial modes are separated according to the different 
polarizations: \H) ^ \H)t\T) and \V) \V) r\R) where 
|r) and \R) refers to the transmitted and reflected spatial 
modes, respectively. The polarization in the reflected 
mode is rotated from \V) to \H), making it identical to 
the transmitted mode. Each mode undergoes the desired 
unitary operations: II(7r/16) for |T) (thus |T^)_r2) and 
H(7r/8) for \R) (thus \V)ri). 

The two spatial modes are arranged in a self- 
stable, displaced Sagnac interferometer configuration 
and mixed into the same PBSe (see Fig. 2). To 
make sure that the measurement projection on |1) of 
the ancilla qubit, the polarization of the transmitted 
mode is rotated 90° (whereas the reflected mode is 
fine as it will naturally lead to a \V) state projec- 
tion). After recombination at the PBSg, the output 
state is: (^ + iy)\+)b\H)Bj\V)R2{sm{Tr/8)\H)A) + (m - 
iy)\~)b\V)Ri\H)R2ismi^/4)\V)A)- 

At the output of the interferometer, a polarizer ori- 
ented at 45° is used to erase the which-path informa- 
tion, which leads to: {{fi + i'){sm{Tr/8)\+)b\H) ri\V) r2) + 
{p - iy)ism{T^/A)\-)b\V)Ri\H)R2))\+)A. HWPs are used 
for the polarization rotations. Small-angle prisms are 
inserted into the interferometer to compensate for un- 
wanted phase shifts. 
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FIG. SI: Observation of Hong-Ou-Mandel interference and four-photon entanglement. A and B show the data of four-photon 
entanglement in the middle stage of running the algorithm (after the PBS3, PBS4 and PBS5). A four-photon Greenberger- 
Horne-Zeilinger (GHZ) state is formed in the state {\H)i\H)2\H)3\H)4 + \V)i\V)2\V)3\V)4)/V2- The multiphoton detection 
probabilities in the H/V measurement basis are listed in A. From left to the right, the x axis is the polarization combinations 
of HHHH, HVHH, HHHV, HVHV, HHVH, HVVH, HHVV, HVVV, VHHH, VVHH, VHHV, VVHV, VHVH, VVVH, VHVV, 
WW. B presents the measured expectation values of the observables X®'', F®*, {X + y)®" and {X - V)®*. The error bars 
denote one standard deviation, deduced from propagated Poissonian counting statistics of the raw detection events. C. Four- 
photon interference as a function of the temporal delay. The red line is a visual guide. As the temporal delay approaches zero, 
the terms \H)i\H)2\H)3\H)4 and V)2! V)3j V)4 gradually become indistinguishable. Thus, the |-|-)i|-|-)2|+)3|+)4 events 
show a maximal enhancement while the |-f)i|+)2|+)3|— )4 events show a dip. The raw visibility of the interference fringe is 
measured to be 0.61(6). Its deviation from the unity is mainly caused by contributions from double pair emissions of parametric 
down conversion. 
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FIG. S2: Reconstructed density matrices of output states |a;). 
Real parts of experimentally obtained (left column) and ideal 
(right column) density matrices are shown for different inputs 
A |6i) = {\H) + \V))/V2, B \b2) = {\H) - \V))/V2, and 
C I63) = \H). Numbers and 1 in a; and y axes represent the 
computational basis of the single qubit. 



